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Abstract 

We propose a new deterministic numerical scheme, based on the discontinuous 
Galerkin method, for solving the Boltzamnn equation for rarefied gases. The new 
scheme guarantees the conservation of the mass, momentum and energy. We avoid 
any stochastic procedures in the treatment of the collision operator of the Boltz- 
mamn equation. 
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1 Introduction 

The classical Boltzmann kinetic equation describes neutral particle transport phenomena. 
Today numerical solutions of the Boltzmann equation are requested to solve problems in 
different fields of real- world applications. There are two classes of computational methods, 
which are used to solve the kinetic equation. In the first of techniques, the well-known 
Direct Simulation Monte Carlo (DSMC) method, the molecular collisions are considered 
on a probabilistic rather than a deterministic basis. The literature of the applications 
of this method is very vast. In the second class, deterministic methods, the Boltzmann 
equation is discretized using a variety of methods and then solved directly or iteratively. 
The computational complexity is high due to the large number of independent variables. 
This heavy computational cost explains why kinetic equations are traditionally simulated 
by the Direct Simulation Monte Carlo methods. As examples of papers dealing with 
deterministic schemes for the Boltzmann equation, we cite some references [T], [TJ, [S] and 
the book [SJ, but there are many other interesting works. 
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In recent years, deterministic solvers to the system, given by Boltzmann equation 
coupled to Poisson equation, describing of electron flow in semiconductors were considered 
in the literature (see, for example, [2], [3], [1]). These methods provide accurate results 
which, in general, agree well with those obtained from DSMC simulations, sometimes 
at a comparable or even less computational time. The discontinuous Galerkin (DG) 
method, which is a finite element method using discontinuous piecewise polynomials as 
basis functions and numerical fluxes based on upwinding for stability, seems to be good for 
solving also kinetic equations. The method has the advantage of flexibility for arbitrary 
unstructured meshes, with a compact stencil, and with the ability to easily accommodate 
arbitrary /jp-adaptivity. For more details about DG scheme for convection dominated 
problems, we refer to the review paper [5]. 

The starting point of this paper is a weak formulation of the Boltzmann equation [U] , 

mo 

If we multiply both sides of the equation by a test function 0(x, £) and we integrate with 
respect to the coordinates x and the velocity then we obtain the equation 

Q -\ f f(t,x,£)4fat)dKd£+ f /Y§£(t,x,fl0(x,fldxd$ 
at Jr3 Jx Jr3 J x tfx 

g(/,/)(t,x,^)0(x^)dxd€, 

X 

where X denotes the x-domain. It is necessary an integration by parts to move the 
derivative with respect to x from the function / to (p. This requires some information on 
the domain X and the test function. We do not describe this step, because, in this paper, 
we will study only the r.h.s of the Boltzmann equation. The test function cj) belongs to a 
suitable chosen finite dimensional space, which is also used to find an approximation of 
the unknown /. 

The plan of the paper is the following. In Section 2, we will introduce the Boltzmann 
equation. Section 2, 3 and 4 will be devoted to the weak form of the collision operator and 
a new modified version. Section 5 will show a simple approximate distribution function 
/ to be used in the framework of the discontinuous Galerkin method. Conclusions and 
future work are given in Section 6. 



2 The Boltzmann equation 

The purpose of this section is to briefly introduce the classical nonlinear Boltzmann equa- 
tion for monatomic gases, to recall well-known properties and to derive simple results. 
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According to the standard notation, the Boltzmamn equation takes the form 



% + ^7T= f I I Wfot.\£',e i )(fZ-ff.)d£.d£'d&. (1) 



The one-particle distribution function / depends on time t, position x and velocity £. 
The kernel W of the collision operator is defined by 

W(t, O = K(n ■ V, | V|) 6(£ O 5(\£ | 2 + |^| 2 - |£f - |£| 2 ) (2) 
where 

n = jfrf[ and V = ^"^- ( 3 ) 

The function K is related to the interaction law between colliding particles. The Dirac 
distributions guarantee momentum and energy conservation during the binary collisions. 
It is immediate to verify that the following symmetry properties 

w{t,t.\e,&) = w(e,&\t,&) and w(t,i.\e,&) = w{t.,£\&e) 

hold. The collision operator of the Boltzmann equation ([[]) is usually written in different 
way, since Dirac distributions are employed and a five-fold integral is derived. In this 
paper we find useful the form of collision term given in ([TJ or the partially reduced 
integral operator, where only the Dirac distribution describing momentum conservation 
is used to reduce the order of integration. 

The integration with respect to the variables and may be performed ab initio in 
the lost term of the collision operator, since the unknown / is not involved. To this aim, 
we must consider the total cross section 

/ / W{t t £.\?,C (4) 

Jk 3 Jk.3 

It is related to the collision frequency. It is easy to prove (see Appendix A) that the 
integral flU) is a function of £ and only through |V| and it is equal to 



\)d», (5) 



where 



k c ((,\v\) = k(^Ji\v\2-(,\v\^ 



3 



2.1 A modified kernel 



The numerical treatment of a kinetic equation by means of finite differences or elements 
requires a bounded domain for the velocity space. We introduce a suitable characteristic 
function in the kernel of the collision operator, such that there exists a bounded domain 
Q so that if /(0,x, £) = for every £ (jL Q and x G X, then /(£, x, £) = for every 
^ fl, x G I and for all time t. Let £ be a positive real number. We define the function 
X £ :l 3 xl 3 ^las follows 

] otherwise 

and 

£) = Xe(t £*) ^1^, £) • (6) 

It is immediate to see that 

/ / W £ (&£.\e,&)ded& = xe(£,£*)S(\V\). (7) 

It is evident that the new modified kernel Ws satisfies the same properties of symmetry 
of true kernel W. Moreover, Ws guarantees that if the particles, before the impact, 
have velocities such that |£| 2 + |£*| 2 < £, then the velocities after the impact satisfy 
the inequality |£'| 2 + |£'| 2 < £ . Using this kernel, for instance, we can choose VL = 
{£ G M 3 : |£| 2 < £}. We can use a smooth function instead of the characteristic function 
Xe to modify the kernel and we obtain the same previous conclusions, but this is useless 
in our numerical approach. 



2.2 Weak form of the collisional operator 

Let (f> : M 3 — > H. be a measurable function. Assuming the existence of the integrals, we 
can recover the well-known result 



w £ (t , o iff* - //*) dL & dc 



R3 JR3 



(8) 



i 

~ 2 
We define 



/ t w(tZ4Z',Oxe(Z,Z*)[cf>' + <f>'*-ci>-<f>*]ff*dZ'dz:d&dZ. 

Jr 3 Jr 3 



w(t & \?, O [<K?) + <K&)\ d? dC 



(9) 
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It is evident that G(<f>; £*) = G((f>; £*, £). We now write the integral OH]) in a meaningful 
form 



f \f f f WeitUZ'iOiffi-ffJd&dZ'dZi 

Jm. 3 Urn. 3 Jr 3 Jr 3 



(10) 



G(<f>; t &) - S(\ V|) [0(£) + 0(C)] Xe (£, &) /(€) /(&) ^* # • (11) 



where, as in the following, to simplify the notation, we omit to write the variables t and 
x, explicitly. In view of numerical calculations of ( ITT]) , we note that the function 



G{<f>;t,&)-S{\V\)[<f>{t) + 



(12) 



which depends on variables and also on the function 0, plays the role of a kernel 

of the integral operator f llip . So, we are aWe to find a reasonable approximation of the 
function ( TJJ| ) /or any /ixed function 0, then we can solve the six-fold integral ( f77]j instead 
of the twelve-fold integral HDi) . There is a simple, but important, case where the function 
()12p is known explicitly. In fact, denoting by ip one of the collision invariants 1, £, |£| 2 , 
we have 

G^;t^)-S(\V\)[m + ^*)] = (13) 



3 The operator G 

It is obvious that simple and explicit expressions of ( I12p are unrealistic in the whole space 
M 3 x M. 3 and for generic 0, apart from the case ( 1T3"]) . The scenario changes drastically, if 
the aim is to find simple approximations of (fT2]) in small domains of the (£, £*) space and 
for a finite set of 0. Since the function 5 is usually given and does not require further 
studies, we must consider only the operator G. 

Let D be a compact set of M 3 such that e M 3 : |£| 2 < £} C £). We look for solutions 
of the Boltzmann equation vanishing for velocities outside the set D. We consider N 
measurable sets C a (a = 1, 2, , .., N) such that 

N 

C a CD Va, C a n^ = (J C a = D . 

01=1 

We denote by Xa the characteristic function on the set C a . 

We are interesting to study G("0 X 7 ; £*)• We need to find a relationship similar to 
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Eq. ( IT21) . Here we have 

N N 



- _ | 7=1 -M 3 JR 3 

JC-, JR. 3 JR 3 JC-y 



7 

N 



£ 

7=1 



JD jr 3 jr 3 jd 



Moreover, 

N 



7=1 

W dZiWe(Z,Z*\Z',Oi>(Z')+ f d£' f d&W e {£,&\e,&)rl>(0 
Jr 3 Jr 3 jd 

d? f dt:We(Z,UZ',£)i>(Z')+ I d? f d£We{t t £.\e,&)il>(&) 
s Jr 3 Jr 3 Jr 3 

[ [ w £ (t^',oma+^(0}d^dc 

Jr 3 Jr 3 

JR 3 JR 3 



Thus 

N 



xe& £*) E G ^ ^ *> = M*) + xs(*. &) 5(| v|) . (14) 

7=1 

3.1 Some formulas 

Taking into account the results of Appendix A, it is easy to see that G(4>; £, £*) can be 
written, for a generic test function 0, as follows 

i / K C (C ■ V, |V|) [0 (§(£ + &) + C) + (§(£ + &) - C)] 5 (|C| 2 - ±|V| 2 ) dC . 
Hence G(ipX-y', £i £*) is the sum of the two integrals, corresponding to the sign + and — , 

\ \ k c (c-v, |v|) v (|(l + &) ± c) x 7 (|(l + £0 ± c) 5 (|C| 2 - ±|v| 2 ) dc . 
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We now consider this operator for the five cases corresponding to the five collision invari- 
ants ip. We obtain the following expressions. 
Case: ip = 1 

\ [ K C (C ■ V, |V|) X y + £*) ± C) S (|C| 2 - ±|V| 2 ) dC . 

Case: ip = $, 

I [ K C {C • V, |V|) + £0 ± C] x 7 m + &) ± C) 5 (|C| 2 - ±|V| 2 ) rfC 

= \ (t + £*) / ^c(c • v, |v|) x 7 + £*) ± c) 5 (|C| 2 - i|v| 2 ) rfc 

±i /* K c (C-V,|V|)Cx 7 (K^ + e)±C)5(|C| 2 -i|V| 2 )rfC. 
Case: ^ = III 2 

Taking into account the Dirac distribution and replacing |C| 2 with ||V| 2 , we have 

[|(€ + £0±c] 2 = |[|€l 2 + l€.| 2 ]±(€ + €.)-c. 

Then 

\ [ K C (C-V,\V\) [i(| + e)±C] 2 x 7 (|(| + e)±C)5(|C| 2 -i|V| 2 )rfC 
= 1 j[\Z\ 2 + lei 2 ] / K C (C ■ V, |V|) x 7 (HZ + &) ± C) 5 (|C| 2 - i|V| 2 ) rfC 

±§(€ + &)- / ^c(C-V,|V|)C X7 (|(| + ^)±c)5(|C| 2 -||V| 2 )rfC. 

In these cases, it is clear that the operator G can be written using the two couples of 
functions 

A±(tH*) = l [ Kc(C-V,|V|) X7 (i(| + e)±C)5(|C| 2 -i|V| 2 )rfC (15) 

1 JR3 

= \ I ^c(CV,|V|)Cx 7 (|(| + e)±C)5(|C| 2 -i|V| 2 )dC. (16) 

1 JR3 

Set 

Mt,**) = + A^(€,€.) and B 7 (£,£*) = B+(£,&) - B;(£,&) , (17) 

for V = 1, I and |£| 2 , G(ip x~n £, £*) writes 
(V> = 1) ^(£,60 

(^=ll| 2 ) -»> ^0^| 2 + I^| 2 ]A 7 (66) + B 7 (66)-(I + ^)- 
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3.2 Approximation of G 



We look for an approximation of G(ip x 7 ; fa for 7 = 1, 2, .., N and any if). In order to 
make clear the problem, we write part of Eq. (TTTj) 



The objective is to find an approximation of the kernel G such that the error introduced 
by the new integral operator is small for a reasonable set of distribution function /. 

We note that G(if)x-y]£,fa is given in terms of the Ay and B 7 (£,£#); so, it is 

sufficient to consider these functions. Let a, (3 and 7 be three positive integers belonging 
to the interval [l,iV]. Suppose M(£) : I 3 4 R a positive function, which represents a 
good candidate of the set of the functions /. We propose these simple approximations 



Xe(t fa Ay^, fa M{$ M(fa « 2 $ 7a/3 X£ (£, S(|V|)M(£) M (&) (19) 
fa B 7 (£, &) M(0 Af (&) « e 7Q ^ £*) I V| 5(| V|) M(£) M(&) (20) 



V£ G C a and ^ t G C^. In Eqs. f|T9|) -f l20|) . $ 7Q ,/3 and the array 7a/ 3 are constant parameters 
to be determined. 

The first step is to consider Eqs. ffT^ - fTSU]) only for a finite number of points £ and 
To this scope, we choose a finite set of well-distributed points in each cell C a . We denote 
these points by £ a ;. 

If Xe{fa fa) fai ~ fa\ S (\€ai ~ fa\) = for every i and j, then we define $ 7Q/3 = and 
7Cr /3 = 0. Otherwise, we use the standard least square method to find the parameters. 
Set y (£,£*) = Xe(€>fa -^(£*)> i n our case the problem is formulated as follows 
find the minimum of the function of $ 7Q/ g 



(18) 



Jr 3 Jm 3 



{ Y (^i, fa) IM^i, fa) - 2 % a p S(\fa - fa\)}} 2 



and the minimum of the function of 7Q/ g 



{Y{fa,fa) \B>y{£ai,fa) - © 7 a/3 \fa ~ fa\ S(\fa - fa\)]} ■ 



The solution is obtained easily and it is given by 



E lY(Z«i,fa)] 2 Mfa,fa) s(\fa - fa\) 



2j2i Y (fa^fa) s (\fa-fa\)] 2 



(21) 
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and 



& 



7a/3 



i,3 



(22) 



8j\ 



i -J 



It is simple matter to show that 



AT N 

^2 *7«/3 = 1 and ^2 ®7"/3 = • 

7=1 



(23) 



7=1 



The proof requires the use of the expressions of G(lx 7 ; £*) and G(£x~/i£,£*) i n terms 

of ^4 7 (^,^*) and B 7 (£, and Eq. (JHJ). It is also evident that the symmetry property 
$70/3 = $7/3a and 7a/ 3 = @ 7/ 3 a hold. 
If we define 

{0 for ifi = 1 

e 7Q/3 for^ = £ (£ e (7 a and & G (7,0 (24) 

©7^ + for ^ = III 2 

then, for every 7 G [l,iV] and V£ G C a , G (7g, we can write the approximation of 
^(^Xyj £) £*) m the following compact way 



x 7 ; 1 £*) xc(€, &) « xc(€, 5(| v|) 



(25) 



4 Approximate weak form of the collisional operator 

Let us consider Eq. ( ITTj) . We have 



f \f f f We{t,£*\e,&)(f f.~ ff*) 



Gty x 7 ; t &) - 5(| v|) [^(0 *,(*) + x 7 (€.)l 



E f dC [ d^[G^ X ^t^)-S(\y\)[mx^)+^)x^)] 

a,/3 a P 
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[ d ^[ G ^x,;t^)-S(\v\)[ij^)5 ja + ^)5^} 



Xs($,$*) /(£) /(£*) ( we now use the approximation (1251) of G) 



ct,/3 a P 

- Mt) &ya + &rp]] xe{£, £0 /(£) /(&) 

5 - *rJ / ^ / d£*S{\V\)m)xe{t(i*)f{€)f{t*) 

+ ^ Z) - / d6.5(|V|)^(€.)xc(€,6)/(€)/(e.) 

a ^3 ^ Ca ./ 

+ ^ 5(|V|) il^^; £,€,)xff(€, 

i J2 [%ap - 5 ia ] f dtf de.S(|V|M£)x5(& &)/(*)/(&) 

+ ~X>^-<U/ ^/ #*s(|V|)^(eOxe& €.)/(*)/(&) 

+ d ^S(\V\)R ia ^;t^)Xe(t^)f^)f(^)- 

n a J C a JCr 



Therefore we have 



Jr 3 Ivr 3 Jr 3 7r 3 

a B C a ^C/3 



This is the proposed model to approximate 

Q(f-f) x 7 (£) HZ) d Z ■ 



Remark 1. If we make the sum with respect to all values of 7, then the r.h.s of Eq. (1261) 
vanishes since ip is a collision invariant. Now, due to the properties f )23|) . also the l.h.s 
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vanishes. Thus the modified collision operator guarantees the conservation of the mass, 
momentum and energy as the true one. 

Remark 2. We have not introduced any constrain on the size and shape of the cell C a ; 
so, there is a great arbitrariness in the decomposition of the domain D. 
Remark 3. We need the constant parameters $ 7a/ 3 and 7Q( g, which depend only on the 
decomposition of D and the function M(£); hence, they be evaluated once at beginning. 
The formulas ( 121 j) and (j22p require to solve the integrals (lTo*j) and (I16p . There is an 
high complexity, due to the presence both the characteristic function x 7 and the Dirac 
distribution. It can be eliminated by means of a transformation of coordinates (from 
cartesian to spherical), by this changes the shape of domain of integration. A very regular 
cell C 7 , having a simple geometry, is usually transformed in a warped domain. In our 
scheme we do not need accurate quadrature formulas, because we have Eqs. (j2"3"|) . These 
can used as a corrector. For instance, we can find the parameters $ 7Q/ 3, use (j2"3"|b for 
checking the goodness and define the new values by means of the simple formula 

[%^) new = — V a, p, 7 . 

A=l 

We can use (1231) 9 to adjust the other parameters 7Q( g. 



5 Approximate distribution functions 

In the framework of the discontinuous Galerkin method, a consistent approximation is 
given using in each cell the same set of functions used as test functions. We denote by 
{ipj : j = 0,1.., 4} the ordered set |£| 2 }. We now assume 

N 

f(t, x, « P(£) \Ba(!> x ) ' Xo(0 • (27) 

We have introduced a weight positive function P(£), which we can choose taking into 
account a reasonable or expected shape of the solution /. The components of the five 
dimensional array r] a (£) are functions, denoted by which are linear combination 

of the collision invariants and such that 

/ p(*)w*)^(e)de = fy (28) 

for every i and j and for each cell C a . The vector functions g a (t, x) (a = 1, 2, .., AT) are 
the new unknowns instead of the distribution function /. It is immediate verify that 

/ /(t,x,£)^(£)d£«<fcj(*,x), (29) 

J Ca. 
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where g a ,j(t, x) are the components of g a (t, x). Using the approximation fl2~T|) . the integrals 
in the r.h.s. of Eq. f[2T)j) become second degree polynomials in the new unknowns g a (£, x). 
The numerical coefficients are constant and determined by solving the integrals 

f d£ f d&S(\V\)mXe(tt*)P(£)P(Z*)r}*M)VpAZ*) (30) 

J Ca. J Cft 



and 

I d£ f d^S(\V\)R, a ^;t^)Xe(t^)P^)P(^)VaA^V^)- (31) 

Also these parameters can be find at the beginning because they depend only on the 
domain decomposition and the weight function P(£). 



6 Conclusions and future work 

We have proposed a numerical model for a deterministic treatment of the collision operator 
in the framework of the discontinuous Galerkin method. The free streaming term of the 
Boltzmann equation was studied [TJ. A numerical application of the proposed technique 
requires the computation of a lot of integrals. This work is required only at beginning. 
We made some very preliminary tests, where we have used is a simple quadrature method. 
For instance, suppose that we must integrate the product of two functions p(z) q(z) where 
z G M m . Let R be the domain of integration and Ri (i = 1,2, ...) a partition of R. We 
denote by suitable point in R^. Then 

p(z)q(z) dz = p(z)q(z)dz^^p(zi) J q(z) dz . 

J R j J Ri j J Ri 

This is useful, whether the integrals of the function q are solved in Ri analytically. The 
function p is the "bad" part of the integral. In our case the function p contains charac- 
teristic functions and *S^(| V|) eventually, and p is a polynomial in the velocity variables. 

The validity of the numerical scheme will be considered in a future work. Some very 
preliminary numerical experiments seem promising. 



7 Appendix A 

We derive a few formulas, which are present in the literature, both to introduce the 
notation and to make the paper self-consistent. Using the properties of Dirac distributions, 
we show how to reduce the order of the integral 

/ / F($,z*,£',&s(z+£*-£'-om 2 +\u 2 -\z?-\z:\ 2 )dz'd£:. (32) 

Jr 3 Jr 3 
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Here F is a generic function so that the following operations are allowed. Introducing the 
change of variables £' = £ — A, = + B, the integral (j3"21) becomes 

/ / £*, I - A, C* + B) 5(A - B) 5(\£\ 2 + |e| 2 - II - A| 2 - |& + B| 2 ) rfA rfS 

and, after a simple integration, 

J / F(£,&,S-A,& + A)5(\A\ 2 -A-V)dA. 

Here we have used the vector V, which is defined in fl3]). Now, performing the further 
change of variable C = — A, it is easy to verify that the integral (|32p is 

\ [ F{t,£.,l{t + £.) + C,±{t + £.)-C)6(\C\ 2 -l\V\ a ) dC. (33) 
One last reduction is available using the Dirac distribution in (13"3"]) . 

7.1 The total cross section 

A simple application of f )33|) allows to reduce the integral 

Jw. 3 Jm.3 

to 

i / ir(V-n c ,|V|)5(|C| 2 -i|V| 2 )dC, (34) 

where 

nc i£-i(£+eo-cr 

With easy calculations and replacing |C| 2 with ||V| 2 , we obtain 

V n c = ^|V| 2 - V- C. 

Therefore we define 

K C (V ■ C, |V|) = ^^i|V|2-V-C, |V|) , 

and we introduce a reference frame in M. 3 and spherical coordinates such that 
C = r(sin ip cos 9, sin ip sin 0, cos ip) and C ■ V = r |V| cos<^ . 
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Hence the integral (J3U) reduces to 



1 /»+0© /*7T /•27T 

-/ dr dip deK c (r\V\cosip, |V|) <5 (r 2 - ±|V| 2 ) - 

2 Jo io io 



7T 



V| / K c (±|V| 2 cos(/?, |V|)shi(^ = ^-|V| / # C (§|V|V, |V|) dp 



7T , 



r sin </? 



1 117-12, 



-1 
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